library(rmarkdown)
library(dplyr)
library(tidyverse) # pour %>%
library(kableExtra) # kbl
library(knitr) # kable
library(ggplot2)
library(plotly) # graphes dynamiques
library(MASS)
library(corrplot) # corrplot
library(GGally) # ggpairs
library(gridExtra) # co-plot pour ggplot
library(ggfortify) # autoplot
library(car) # Tests type II
library(rsample) # create training/test sets
library(corrplot)Corrélation : La corrélation mesure la relation statistique entre deux variables, indiquant si et dans quelle mesure elles varient ensemble.
Causalité : La causalité désigne la relation de cause à effet où un événement (cause) entraîne un autre événement (effet), établissant une relation de dépendance directe.
⚠️ La corrélation ne prouve pas la causalité.
Ici on observe une corrélation entre deux variables “Plus d’éducation = moins de chômage”, mais cela ne signifie pas nécessairement qu’il y a une relation de cause à effet (causalité) entre ces variables :
Corrélation Spurieuse : Des facteurs non mesurés tels que l’économie, la politique ou le social peuvent influencer à la fois l’éducation et le chômage sans qu’il y ait un lien direct entre eux.
Direction de la Causalité Inverse : Il est possible que ce soit une économie prospère qui encourage davantage l’éducation plutôt que l’éducation seule qui réduise le chômage.
Facteurs Confondants : Des variables telles que la situation économique globale ou les politiques gouvernementales peuvent influencer à la fois l’éducation et le chômage, créant ainsi une corrélation apparente.
Causalité Complexes : La relation entre l’éducation et le chômage peut être complexe, impliquant des interactions multiples et non linéaires avec d’autres facteurs, ce qui rend difficile l’attribution d’une causalité basée uniquement sur la corrélation observée.
En résumé, bien qu’il puisse exister une corrélation entre le niveau d’éducation et le taux de chômage, d’autres facteurs doivent être pris en compte pour déterminer s’il existe réellement une relation de cause à effet entre ces variables. Des études plus approfondies, des analyses de causalité et la prise en compte de variables confondantes sont nécessaires pour établir des liens de causalité robustes.
Dans ce cas on observe une corrélation linéaire entre deux variables telles que le nombre de prix Nobel dans un pays et la consommation de chocolat, il est peu probable que cela implique une relation de causalité directe entre ces variables. Cette situation illustre le principe selon lequel la corrélation ne prouve pas la causalité.
En résumé, bien que la corrélation puisse être observée entre le nombre de prix Nobel dans un pays et la consommation de chocolat, il est important de ne pas conclure à une relation de causalité directe sur la seule base de la corrélation. Des études plus approfondies, une analyse des mécanismes causaux potentiels et la prise en compte d’autres variables sont nécessaires pour évaluer la nature de la relation entre ces variables.
Ce jeu de données (Air Breizh, Summer 2001) contient 112 observations et 13 variables
Corrigeons la déclaration erronnée des variables:
#library(dplyr)
ozone1$Date <- as.Date(as.factor(ozone1$Date),format = "%Y%m%d")
ozone1$vent<-as.factor(ozone1$vent)
ozone1$pluie<-as.factor(ozone1$pluie)
ozone1 <- ozone1 %>% mutate(across(where(is.character), as.numeric))
ozone1 <- ozone1 %>% mutate(across(where(is.integer), as.numeric))
paged_table(ozone1)💡 Bon à savoir !
Pour changer la nature en date, il faut d’abord transformer en facteur, puis préciser le format.
Ici, across est utilisé avec
where(is.character) pour cibler uniquement les colonnes qui
sont de type caractère (chr) et appliquer as.numeric à ces
colonnes.
But Prédire la concentration maximum d’ozone dans
l’air (maxO3). Interessons nous dans un premier temps
uniquement à la Températures relevées à 12h (T12). Trouver
une fonction \(f\) telle que
maxO3\(\approx
f\)(T12)
avec
T12 Températures
relevées à 12h.maxO3 :
Concentration maximum d’ozone dans l’air.Comment choisir \(f\)
\(\Rightarrow\) Trouver la fonction qui minimise la somme des carrés des distances de la fonction aux points (\(T12_i,maxO3_i\)).
Intéressons-nous au modèle linéaire, c’est-à-dire que
maxO3\(\approx\)
fonction affine de T12 \(\Leftrightarrow\) \(Y=f(X)+\epsilon=\beta_0+\beta_1X+\epsilon\)
T12 : variable
réponse, variable à expliquer,..maxO3 : variable
explicative, régresseur, prédicteur, covariable,..Problème \((\beta_O,\beta_1)\) sont inconnus!
Solution les estimer à partir de nos données, de notre échantillon de \(n\) observations \((X_i,Y_i)_{i=1,\cdots,n}\).
Modèle linéaire simple:
\[Y_i=f(X_i)+\epsilon_i=\beta_0+\beta_1X_i+\epsilon_i,\qquad i=1,\cdots,n\]
[Hypothèses:]
On peut résumer ces hypothèses par \(\epsilon_i\overset{iid}{\sim}\mathscr N(0,\sigma^2)\)
[Prédiction le modèle estimé:] \[\widehat Y_i=\widehat f\big(X_i\big)= \widehat\beta_0+\widehat\beta_1X_{i}\]
Pour trouver \((\widehat\beta_0,\widehat\beta_1)\), on minimise la somme des erreurs aux carrés à partir du modèle \[Y_i=f(X_i)+\epsilon_i=\beta_0+\beta_1X_i+\epsilon_i,\qquad i=1,\cdots,n\] C’est à dire \((\widehat\beta_0,\widehat\beta_1)\) sont les arguments qui minimise la somme des carrés résiduels \[ (\widehat\beta_0,\widehat\beta_1)=\underset{(\beta_0,\beta_1)}{\arg\min}\sum_{i=1}^n(Y_i- f(X_i))^2=\underset{(\beta_0,\beta_1)}{\arg\min}\sum_{i=1}^n(Y_i-\beta_0-\beta_1X_i)^2=\underset{(\beta_0,\beta_1)}{\arg\min}\sum_{i=1}^n\epsilon_i^2 \] On appelle cet estimateur, l’estimateur des moindres carrés, qui est également dans ce cadre l’estimateur du maximum de vraisemblance. Avec R:
| Estimation | |
|---|---|
| (Intercept) | -27.419636 |
| T12 | 5.468685 |
Traçons cette droite, appellée la droite des moindres
carrées.
ggplot(ozone1)+ aes(x = T12, y =maxO3)+ geom_point(col='red', size = 0.5)+ geom_smooth(method = "lm",se = FALSE)Si maintenant on considère nos observations en prenant en compte le
facteur ventet que l’on souhaite faire régression
différente pour chaque modalité.
p1<-ggplot(ozone1)+ aes(x = T12, y =maxO3)+geom_point(size = 0.5)
p2<-p1+ geom_smooth(method = "lm",se = FALSE)
p3<-ggplot(ozone1)+ aes(x = T12, y =maxO3, color=vent)+geom_point(size = 0.5)
p4<- p3+guides(color = FALSE)
p5<-p3+ geom_smooth(method = "lm",se = FALSE)
subplot(p1,p2,p4,p5,nrows = 2)Alors il faut définir autant de modèle linéaire du type
\[Y_{i}=f(X_i)+\epsilon_i=\beta_{0}+\beta_{1}X_{i}+\epsilon_{i}, i=1,\cdots,n\] qu’il y a de modalités \(k\).
Modèle linéaire (Forme régulière):
\[Y_{ik}=f(X_{ik})+\epsilon_{ik}=\beta_{0k}+\beta_{1k}X_{ik}+\epsilon_{ik},\qquad i=1,\cdots,n_k, \ k=1,\dots,4\] où
vent
.vent
.Mais on pourrait imaginer qu’il existe un intercept et une pente de reférence identiques à tous
vent.T12et le facteur
vent.Modèle linéaire (Forme singulière):
Pour tout \(i=1,\cdots,n_k,\) et tout \(k=1,\dots,4\) \[Y_{ik}=f(X_{ik})+\epsilon_{ik}=\mu+\alpha_k+(b+c_k)X_{ik}+\epsilon_{ik}\]
⚠️ Ce modèle n’est pas identifiable,il existe une infinité de solutions. Pour le rendre identifiable, il faut ajouter des contraintes.
Par défaut R ajoute la contrainte \(\alpha_1=c_1=0\) pour “exhiber” une solution.
lm(maxO3~T12+vent+T12:vent, data=ozone1)
ou bien
lm(maxO3~T12*vent, data=ozone1)
mod_defaut<-lm(maxO3~T12*vent, data=ozone1)
mod_defaut$coefficients%>%kbl(col.names=c("Estimation"))| Estimation | |
|---|---|
| (Intercept) | 8.735714 |
| T12 | 4.092281 |
| ventNord | -24.902135 |
| ventOuest | -57.366714 |
| ventSud | -43.195199 |
| T12:ventNord | 1.006046 |
| T12:ventOuest | 2.232719 |
| T12:ventSud | 1.680646 |
ventEst
n’apparait pas car R a pris par défault \(\widehat\alpha_1=0\)T12:ventEst
n’apparait pas car R a pris par défault \(\widehat b_1=0\)Retrouver le modèle régulier il suffit sous R d’indiquer que \(\mu=b=0\)
lm(maxO3~-1+T12+vent, data=ozone1)
mod_reg<-lm(maxO3~-1+T12+vent, data=ozone1)
library(kableExtra)
mod_reg$coefficients%>%kbl(col.names=c("Estimation"))| Estimation | |
|---|---|
| T12 | 5.479838 |
| ventEst | -24.107753 |
| ventNord | -23.821255 |
| ventOuest | -30.814974 |
| ventSud | -27.504905 |
Ìntercept n’apparait pas car R a pris
par défault \(\widehat\mu=0\)T12 n’apparait
pas car R a pris par défault \(\widehat b=0\)💡 Bon à savoir !
Quel que soit le modèle et les contraintes choisis, la prédiction sera la même ; seule l’interprétation des paramètres peut différer.
But Prédire la concentration maximum d’ozone dans
l’air (maxO3) en fonction des variables de notre jeu de
données à partir d’un modèle linéaire.
Dateset.seed(2024)
ozone_split <- rsample::initial_split(ozone1bis, prop = 0.8)
ozone<- rsample::training(ozone_split)
ozone_test <- rsample::testing(ozone_split)maxO3 :
Concentration maximum d’ozone dans l’air.T9
Températures relevées à 9h.T12
Températures relevées à 12h.T15
Températures relevées à 15h.Ne9
Variable synthétique nebuleuse à 9h.Ne12
Variable synthétique nebuleuse à 12h.Ne15
Variable synthétique nebuleuse à 15h.Vx9
Variable synthétique vent à 9h.Vx12
Variable synthétique vent à 12h.Vx15
Variable synthétique ventà 15h.maxO3v
.vent
Facteur vent à 4 modalités (Est, Nord,
Ouest, Sud).pluie
Facteur pluie à 2 modalités (Pluie, Sec).maxO3\(=f\)(T12,T9,T12,T15,Ne9,Ne12,Ne15,Vx9,Vx12,Vx15,maxO3v,vent,pluie)\(+\epsilon\)
\[ Y=f\big(X^{(1)},\cdots,X^{(12)}\big)+\epsilon \]
où les mêmes hypothèses sont faites:
[Hypothèses:]
Prédiction \(\widehat Y_i=\widehat f\big(X^{(1)},\cdots,X^{(12)}\big)\)
Si on ne considère aucune interactions alors sous R le modèle par défaut est le modèle singulier et s’écrit
lm(maxO3~T9+T12+T15+Ne9+Ne12+Ne15+Vx9+Vx12+Vx15+maxO3v+vent+pluie,data=ozone)
ou bien plus rapidement
lm(maxO3~.,data=ozone)
Déclarons notre modèle complet mod_completet affichons
la sortie summary(mod_complet)
##
## Call:
## lm(formula = maxO3 ~ ., data = ozone)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.518 -8.506 -0.731 7.450 44.092
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.2487 18.1291 1.227 0.2236
## T9 1.5073 1.4334 1.052 0.2964
## T12 0.6172 1.8575 0.332 0.7406
## T15 0.2322 1.4355 0.162 0.8720
## Ne9 -2.6971 1.1829 -2.280 0.0255 *
## Ne12 0.4995 1.7572 0.284 0.7770
## Ne15 -0.6239 1.2956 -0.482 0.6316
## Vx9 0.4761 1.0939 0.435 0.6647
## Vx12 0.6873 1.4162 0.485 0.6289
## Vx15 0.4899 1.0771 0.455 0.6505
## maxO3v 0.3822 0.0775 4.931 4.87e-06 ***
## ventNord -0.8189 7.8379 -0.104 0.9171
## ventOuest 3.1013 9.4555 0.328 0.7438
## ventSud 4.6715 8.4367 0.554 0.5814
## pluieSec 2.9572 4.1185 0.718 0.4750
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.88 on 74 degrees of freedom
## Multiple R-squared: 0.7505, Adjusted R-squared: 0.7033
## F-statistic: 15.9 on 14 and 74 DF, p-value: < 2.2e-16
Regardons de plus près
On observe 4 colonnes
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.2487 18.1291 1.227 0.2236
T9 1.5073 1.4334 1.052 0.2964
T12 0.6172 1.8575 0.332 0.7406
T15 0.2322 1.4355 0.162 0.8720
Ne9 -2.6971 1.1829 -2.280 0.0255 *
Ne12 0.4995 1.7572 0.284 0.7770
Ne15 -0.6239 1.2956 -0.482 0.6316
Vx9 0.4761 1.0939 0.435 0.6647
Vx12 0.6873 1.4162 0.485 0.6289
Vx15 0.4899 1.0771 0.455 0.6505
maxO3v 0.3822 0.0775 4.931 4.87e-06 ***
ventNord -0.8189 7.8379 -0.104 0.9171
ventOuest 3.1013 9.4555 0.328 0.7438
ventSud 4.6715 8.4367 0.554 0.5814
pluieSec 2.9572 4.1185 0.718 0.4750
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Estimate : correspond aux \(\widehat\beta_j\) l’estimation des
paramètres \(\beta_j\) par maximum de
vraissemblance/moindres carrésStd. Error : Racine carré de l’erreur d’estimation des
\(\beta_j\)t value: Valeur de la statistique de test du test :
\(H_0\, : \ \beta_j=0 \ \ vs \ H_1\, : \ \
\beta_j\neq 0\)Pr(>|t|) : \(p\)-value du test (de Student)
précédent.⚠️ Attention !
Ne pas rejeter \(H_0\) \(\large{\neq}\) accepter \(H_1\). Ne signifie pas que la variable n’est pas pertinente (écarter une variable est ici prématuré).
💡 Bon à savoir ! Pour afficher uniquement ces 4 colonnes
Dans la sortie du summary(mod_complet) nous avons
également la ligne suivante :
Test de Fisher global: Teste la pertinence des variables dans leur ensemble (sous \(H_0\) modèle réduit à l’intercept)}.
F-statistic: 15.9 : Valeur de la statistique de test:
\(H_0\, : \ \forall j: \ \beta_j=0 \ \ vs \
H_1\, : \overline{H_0}\)p-value: < 2.2e-16 : \(p\)-value du test. Ici, on rejette
\(H_0\), il existe au moins une
varibale pertinente.Dans la sortie du summary(mod_complet) nous avons
également les lignes suivantes :
Residual standard error: 14.88 on 74 degrees of freedom
Multiple R-squared: 0.7505, Adjusted R-squared: 0.7033 Residual standard error : Estimation \(\widehat\sigma\) du niveau de bruit \(\sigma\) (écart type) défini dans les
hypothèses du modèle. Ici \(\widehat\sigma=14.88\), très élevéMultiple R-squared : coefficient determination \(R^2\)Adjusted R-squared : coefficient determination ajusté
\(R^2_a\)Interprétation des coeff de determination: Ils sont
compris entre 0 et 1. La valeur 0 indique que le modèle ne
colle pas aux données, la valeur 1indique que le modèle
parfaitement aux données.
⚠️ Attention ! Plus le nombre de variables explicatives est grand, plus le \(R^2\) est grand! Le coefficient \(R^2_a\) est moins impacté par ce nombre (=la dimension)
Pour finir, dans la sortie du summary(mod_complet) nous
avons également les lignes suivantes :
Residuals : Les résidus estimés (les théoriques \(\epsilon_i\) sont inconnus) définis de la
façon suivantes: \[\widehat\epsilon_i=Y_i-\widehat Y_i=Y_i-\widehat
f\big(X^{(1)},\cdots,X^{(12)}\big)\]📌 Pour valider le modèle (c’est-à-dire les hypothèses), il faudra faire une analyse des résidus. Mais avant de valider le modèle, avons-nous sélectionner le “meilleur” modèle ? Toutes les variables sont-elles pertinentes ?
library(dplyr)
library(corrplot)
ozone_numerique<- ozone %>% select_if(is.numeric)
correlation_matrix<-cor(ozone_numerique)
corrplot(correlation_matrix, order = 'hclust',addrect = 3)Les variables numériques sont très corrélées entre elles!
Méthodes pas-à-pas:
On choisit un critère (\(R^2_a\), AIC (par défaut), BIC, \(C_p\)-Mallows,…) ET une des 3 méthodes pas-à-pas suivantes
Méthode Backward : On part du modèle de taille \(p\) variables explicatives et on enlève 1 à 1 les variables en comparant les modèles 2 à 2 selon un critère.
Méthode Forward : On part du modèle de taille \(1\) réduit à l’intercept et on ajoute 1 à 1 les variables en comparant les modèles 2 à 2 selon un critère.
Méthode Both: mixte de 2 méthodes. On part de l’intercept et on ajoute/enlève les variables 1 à 1 et on compare selon un critère.
Arrêt lorsque le critère n’est plus améliorable .
mod_0=lm(maxO3~1,data=ozone)
mod_forw<-step(mod_0, maxO3~T9+T12+T15+Ne9+Ne12+Ne15+Vx9+Vx12+Vx15+maxO3v+vent+pluie,data=ozone,trace=F,direction=c('forward'))
mod_forw##
## Call:
## lm(formula = maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone)
##
## Coefficients:
## (Intercept) T12 maxO3v Ne9 Vx12
## 20.6041 2.2690 0.3923 -2.5949 1.0888
mod_both<-step(mod_0,maxO3~T9+T12+T15+Ne9+Ne12+Ne15+Vx9+Vx12+Vx15+maxO3v+vent+pluie,data=ozone, trace=F,direction=c("both"))
mod_both##
## Call:
## lm(formula = maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone)
##
## Coefficients:
## (Intercept) T12 maxO3v Ne9 Vx12
## 20.6041 2.2690 0.3923 -2.5949 1.0888
##
## Call:
## lm(formula = maxO3 ~ T9 + Ne9 + Vx9 + maxO3v, data = ozone)
##
## Coefficients:
## (Intercept) T9 Ne9 Vx9 maxO3v
## 24.4931 2.6017 -2.9336 1.6811 0.3845
Les 2 méthodes donne le même modèle, c’est un hasard. Mais comparer au modèle complet et celui restreint à l’intercept, quel est le meilleur model?
Notre modèle réduit à 4 variables + l’intercept est meilleur selon le critère AIC.
##
## Call:
## lm(formula = maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone)
##
## Residuals:
## Min 1Q Median 3Q Max
## -54.128 -9.989 -0.907 8.630 43.605
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.60408 12.63631 1.631 0.10673
## T12 2.26903 0.54497 4.164 7.55e-05 ***
## maxO3v 0.39230 0.06874 5.707 1.68e-07 ***
## Ne9 -2.59485 0.87285 -2.973 0.00385 **
## Vx12 1.08876 0.70288 1.549 0.12514
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.33 on 84 degrees of freedom
## Multiple R-squared: 0.7375, Adjusted R-squared: 0.725
## F-statistic: 59.01 on 4 and 84 DF, p-value: < 2.2e-16
Notons que si le \(R^2\) a diminué le \(R^2_a\) a augmenté. Les erreurs d’estimation ont également diminuées. Le modèle est meilleur.
Test de type I
Ligne à ligne un test est fait avec sous \(H_0\) un modèle avec \(d+1\) variables vs sous \(H_1\) un modèle avec \(d+2\). On ajoute les variables 1 à une, pour tester leur pertinence
intercept vs \(H_1\) intercept+T9intercept+T9 vs \(H_1\) intercept+T9+T12intercept+T9+T12 vs \(H_1\)
intercept+T9+T12+T15intercept+T9+T12+T15 vs \(H_1\)
intercept+T9+T12+T15+Ne9intercept+T9+T12+T15+Ne9 vs \(H_1\)
intercept+T9+T12+T15+Ne9+Ne12👀 Ici on garde les variables :
T9, T12, Ne9, Vx9et
maxO3v. Définissons ce modele
Test de type II Ligne à ligne un test est fait avec sous \(H_0\) un modèle avec \(p\) variables vs sous \(H_1\) un modèle avec \(p+1\). On retire une variable pour tester sa pertinence
T9T12T15Ne9Ne12👀 Ici on garde les variables :
Ne9 et maxO3v. Définissons ce modele
Comparons nos modeles
Le modèle mod_both reste le meilleur suivant le critère
AIC.
Hypothèse [P1]: Erreurs sont centrées/le modèle est linéaire.
📌 Les fitted sont les
réponses prédites \(\widehat Y_i\) par
le modèle.
span class=“solution”>👀 Ici [P1] validée
Hypothèse [P2]: Variance est homoscédastique.
span class=“solution”>👀 Ici [P2] validée
## [1] 0.006507032
span class=“solution”>👀 Ici \(p\)-value<5 % donc [P2] non validée
Hypothèse [P3]: Erreurs sont non corrélées.
## [1] 0.264
span class=“solution”>👀 Ici \(p\)-value>5 % donc [P3] validée
Hypothèse [P4]: Erreurs sont gaussiennes.
Si dans le plot Normal Q-Q:n, les points sont alignés autour de la première bissectrice alors [P4] est validée.
Si il est difficile de valider alors faire test de Shapiro pour décider. \(H_0\): “Gaussien”
## [1] 0.01259202
span class=“solution”>👀 Ici \(p\)-value>5 % donc [P4] non validée
Outliers à fort résidus:
Le modèle peut échoué à prédire correctement \(Y_i\) \(\Rightarrow\) grande différence entre \(|Y_i-\widehat Y_i|=|\epsilon_i|>2\). On peut identifier ces observations sur le Studentized-plot.
2 observations sont des outilers à fort résidus : 31et
37
Outliers à fort effet levier sur eux-même:
Si une observation avec grande influence sur sa propre estimation (\(\widehat Y_i\)), on dit qu’elle à un fort effet levier sur sa propre estimation : \(\Rightarrow\) la hat-value>0.5. On peut identifier ces observations sur le Hat-plot
Aucune observation n’a de fort impact sur sa propre estimation.
Outliers à fort effet levier sur le modèle:
Si une observation a une grande influence sur le modèle lui-même (\(\widehat \beta_j\)), on dit qu’elle à un fort effet levier sur le modèle : \(\Rightarrow\) si sa distance de Cook est grande. On peut identifier ces observations sur le Cook-plot
Aucune observations ne semble avoir une distance très grande.
L’observation 31 a la distance de cook la plus grande.
Mais doit-on pour autant retirer ces points ? Tous les ouliers sont-ils à retirer ?
Le critère de Bonferroni est une méthode d’ajustement pour tenir compte du fait que plusieurs tests sont effectués simultanément. Il contrôle le taux d’erreur de famille, c’est-à-dire la probabilité de commettre au moins une erreur de type I parmi tous les tests effectués.
Dans le contexte de influenceIndexPlot avec
vars = "Bonf", les \(p\)-values associées à chaque
observation sont ajustées selon le critère de Bonferroni. Les
observations qui ont des \(p\)-values ajustées inférieures à
un seuil spécifié sont considérées comme ayant une influence
significative sur le modèle.
Lorsque vous utilisez influenceIndexPlot avec l’argument
vars = "Bonf", cela signifie que vous souhaitez inclure des
lignes dans le graphique des indices d’influence qui correspondent aux
observations ayant une influence significative sur le modèle, selon le
critère de Bonferroni.
Les outliers problématiques ont une \(p\)-value de Bonferroni \(<\) seuil On peut afficher le
Bonferroni-plot
Il est important de noter que l’interprétation des p-values ajustées dépend du seuil alpha choisi. Un seuil alpha plus strict (par exemple, 0.01) conduira à moins d’observations étiquetées comme ayant une influence significative, tandis qu’un seuil alpha plus permissif (par exemple, 0.05) conduira à plus d’observations étiquetées.
En résumé, lorsque vous utilisez vars = "Bonf" avec
influenceIndexPlot, les observations étiquetées ont des
p-values ajustées selon le critère de Bonferroni qui sont inférieures au
seuil alpha spécifié, ce qui indique une influence significative sur le
modèle de régression. Adjuster le seuil alpha permet de contrôler la
sensibilité du critère d’inclusion des observations dans le graphique
des indices d’influence.
Mais le plus simple est d’acceder directement aux \(p\)-value des points problématiques. Pour afficher les \(p\)-\(values\) des points jugés douteux par R, il suffit d’écrire
## rstudent unadjusted p-value Bonferroni p
## 31 -4.418044 2.9857e-05 0.0026573
Dans notre exemple: l’observation 31 a la \(p\)-value au sens de Bonferroni.
mod_complet<-lm(maxO3~., data=ozone)
mod_F_compl<-lm(maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone)
mod_F_ss_outlier<-lm(maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone[-c(31),])Ozone_actual<-ozone_test$maxO3
Results<-as.data.frame(Ozone_actual)
Results$Modele_complet<-predict(mod_complet,newdata=ozone_test)
Results$Modele_Final<-predict(mod_F_compl,newdata=ozone_test)
Results$Modele_Final_ss_outlier<-predict(mod_F_ss_outlier,newdata=ozone_test)
Resultslibrary(ggplot2)
# Plot pour le modèle complet
p1 <- ggplot(Results, aes(x = Ozone_actual, y = Modele_complet)) +
geom_point(color = "#31AA41") +geom_smooth(method = "lm", se = FALSE,color = "#31AA41") +ylab("Modele complet")
# Plot pour le modèle avec outlier
p2 <- ggplot(Results, aes(x = Ozone_actual, y = Modele_Final)) +geom_point(color = "#AA0012") +geom_smooth(color = "#AA0012",method = "lm", se = FALSE) + ylab("Modele avec outlier")
# Plot pour le modèle sans outlier
p3 <- ggplot(Results, aes(x = Ozone_actual, y = Modele_Final_ss_outlier))+
geom_point(color = "#4200AA") +geom_smooth(color = "#4200AA",method = "lm", se = FALSE) +ylab("Modele ss outlier")
# Affichage des plots côte à côte
grid.arrange(p1, p2, p3, nrow = 2) Calculons l’erreur moyenne faite par chacun de ces modèles sur ces données non vues
library(Metrics)
Erreur_mod_complet<-rmse(ozone$maxO3,predict(mod_complet,newdata=ozone))
Erreur<-as.data.frame(Erreur_mod_complet)
Erreur$Modele_Final<-rmse(ozone$maxO3,predict(mod_F_compl,newdata=ozone))
Erreur$Modele_Final_ss_outlier<-rmse(ozone$maxO3,predict(mod_F_ss_outlier,newdata=ozone))
Erreur_mod_complet<-rmse(ozone_test$maxO3,predict(mod_complet,newdata=ozone_test))
Erreur_test<-as.data.frame(Erreur_mod_complet)
Erreur_test$Modele_Final<-rmse(ozone_test$maxO3,predict(mod_F_compl,newdata=ozone_test))
Erreur_test$Modele_Final_ss_outlier<-rmse(ozone_test$maxO3,predict(mod_F_ss_outlier,newdata=ozone_test))
rbind.data.frame(train=Erreur,test=Erreur_test)library(randomForest)
library(Metrics)
# Créer un modèle de Random Forest
mod_RF <- randomForest(maxO3~.,data =ozone)
#c(rmse(ozone$maxO3),predict(mod_RF,newdata=ozone))Erreur_mod_RF<-rmse(ozone$maxO3,predict(mod_RF,newdata=ozone))
Erreur_mod_RF_test<-rmse(ozone_test$maxO3,predict(mod_RF,newdata=ozone_test))
paste0('RF train: ',Erreur_mod_RF,' RF test: ' ,Erreur_mod_RF_test)## [1] "RF train: 6.45839847004582 RF test: 14.6830507786875"
👀 ICI les RF font de l’overfitting, le
rmsesur le train et le test sont trop différent.
📌 Utilisation des RF pour faire de la selection. Utile lorsque le nombre de variables potentiellement explicatives est important.
Si on prend le même nombre de variables que dansq le ML, c’est-à-dire 4, les RF n’améliorent pas notre prédiction
# Obtenir l'importance des caractéristiques
B<-vip(mod_RF,num_features = 4)
B<-B$data$Variable
ozone_RF<- ozone[,c('maxO3',B )]
mod_LM_RF<-lm(maxO3~.,data = ozone_RF)
Erreur_mod_LM_RF<-rmse(ozone$maxO3,predict(mod_LM_RF,newdata=ozone))
Erreur_mod_LM_RF_test<-rmse(ozone_test$maxO3,predict(mod_RF,newdata=ozone_test))
paste0('LM_RF train: ',Erreur_mod_LM_RF,' LM_RF test: ' ,Erreur_mod_LM_RF_test)## [1] "LM_RF train: 15.6072400107394 LM_RF test: 14.6830507786875"
Reprenons pour la suite le jeu de données cars de la
librairie ‘MASS’ qui contient
50 Données et 2 variables : speed:la vitesse en miles
par heure et dist la distance nécessaire à l’arrêt d’un
véhicule.
library(ggplot2)
library(gridExtra)
plot_p<-ggplot(cars, aes(x = speed, y = dist)) + geom_point( size=0.85)
plot1 <- plot_p + geom_smooth( color="#FF9900", method = "lm" )+ labs(title = "Régression linéaire")
plot2 <-plot_p+ geom_smooth(color="#4169E1", method = "lm", formula = y ~ poly(x, 2, raw = TRUE)) +labs(title = "Polynôme de degré 2")
grid.arrange(plot1, plot2, nrow = 1)Définissons ces modèles
Regardons le plot des résidus
La courbe bleue à une forme parabolique. Essayons alors d’exprimer
disten fonction d’un polynôme d’ordre 2
##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.069 -9.525 -2.272 9.215 43.201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.5791 6.7584 -2.601 0.0123 *
## speed 3.9324 0.4155 9.464 1.49e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438
## F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12
##
## Call:
## lm(formula = dist ~ poly(speed, 2), data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.720 -9.184 -3.188 4.628 45.152
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.980 2.146 20.026 < 2e-16 ***
## poly(speed, 2)1 145.552 15.176 9.591 1.21e-12 ***
## poly(speed, 2)2 22.996 15.176 1.515 0.136
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.18 on 47 degrees of freedom
## Multiple R-squared: 0.6673, Adjusted R-squared: 0.6532
## F-statistic: 47.14 on 2 and 47 DF, p-value: 5.852e-12
Regardons le plot des résidus pour valider le modèle
modele_lineaire<-rmse(cars$dist,mod$fitted.values)
Erreur<-as.data.frame(modele_lineaire)
Erreur$modele_polynomial<-rmse(cars$dist,mod2$fitted.values)
rownames(Erreur) <- "Erreur"
Erreur